<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
 <head>
  <title>slaebz.f</title>
 <meta name="generator" content="emacs 21.3.1; htmlfontify 0.20">
<style type="text/css"><!-- 
body { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.default   { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.default a { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
span.comment   { color: rgb(178, 34, 34);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.comment a { color: rgb(178, 34, 34);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
 --></style>

 </head>
  <body>

<pre>
      SUBROUTINE <a name="SLAEBZ.1"></a><a href="slaebz.f.html#SLAEBZ.1">SLAEBZ</a>( IJOB, NITMAX, N, MMAX, MINP, NBMIN, ABSTOL,
     $                   RELTOL, PIVMIN, D, E, E2, NVAL, AB, C, MOUT,
     $                   NAB, WORK, IWORK, INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  -- LAPACK auxiliary routine (version 3.1) --
</span><span class="comment">*</span><span class="comment">     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
</span><span class="comment">*</span><span class="comment">     November 2006
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Scalar Arguments ..
</span>      INTEGER            IJOB, INFO, MINP, MMAX, MOUT, N, NBMIN, NITMAX
      REAL               ABSTOL, PIVMIN, RELTOL
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Array Arguments ..
</span>      INTEGER            IWORK( * ), NAB( MMAX, * ), NVAL( * )
      REAL               AB( MMAX, * ), C( * ), D( * ), E( * ), E2( * ),
     $                   WORK( * )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Purpose
</span><span class="comment">*</span><span class="comment">  =======
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  <a name="SLAEBZ.22"></a><a href="slaebz.f.html#SLAEBZ.1">SLAEBZ</a> contains the iteration loops which compute and use the
</span><span class="comment">*</span><span class="comment">  function N(w), which is the count of eigenvalues of a symmetric
</span><span class="comment">*</span><span class="comment">  tridiagonal matrix T less than or equal to its argument  w.  It
</span><span class="comment">*</span><span class="comment">  performs a choice of two types of loops:
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  IJOB=1, followed by
</span><span class="comment">*</span><span class="comment">  IJOB=2: It takes as input a list of intervals and returns a list of
</span><span class="comment">*</span><span class="comment">          sufficiently small intervals whose union contains the same
</span><span class="comment">*</span><span class="comment">          eigenvalues as the union of the original intervals.
</span><span class="comment">*</span><span class="comment">          The input intervals are (AB(j,1),AB(j,2)], j=1,...,MINP.
</span><span class="comment">*</span><span class="comment">          The output interval (AB(j,1),AB(j,2)] will contain
</span><span class="comment">*</span><span class="comment">          eigenvalues NAB(j,1)+1,...,NAB(j,2), where 1 &lt;= j &lt;= MOUT.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  IJOB=3: It performs a binary search in each input interval
</span><span class="comment">*</span><span class="comment">          (AB(j,1),AB(j,2)] for a point  w(j)  such that
</span><span class="comment">*</span><span class="comment">          N(w(j))=NVAL(j), and uses  C(j)  as the starting point of
</span><span class="comment">*</span><span class="comment">          the search.  If such a w(j) is found, then on output
</span><span class="comment">*</span><span class="comment">          AB(j,1)=AB(j,2)=w.  If no such w(j) is found, then on output
</span><span class="comment">*</span><span class="comment">          (AB(j,1),AB(j,2)] will be a small interval containing the
</span><span class="comment">*</span><span class="comment">          point where N(w) jumps through NVAL(j), unless that point
</span><span class="comment">*</span><span class="comment">          lies outside the initial interval.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Note that the intervals are in all cases half-open intervals,
</span><span class="comment">*</span><span class="comment">  i.e., of the form  (a,b] , which includes  b  but not  a .
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  To avoid underflow, the matrix should be scaled so that its largest
</span><span class="comment">*</span><span class="comment">  element is no greater than  overflow**(1/2) * underflow**(1/4)
</span><span class="comment">*</span><span class="comment">  in absolute value.  To assure the most accurate computation
</span><span class="comment">*</span><span class="comment">  of small eigenvalues, the matrix should be scaled to be
</span><span class="comment">*</span><span class="comment">  not much smaller than that, either.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  See W. Kahan &quot;Accurate Eigenvalues of a Symmetric Tridiagonal
</span><span class="comment">*</span><span class="comment">  Matrix&quot;, Report CS41, Computer Science Dept., Stanford
</span><span class="comment">*</span><span class="comment">  University, July 21, 1966
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Note: the arguments are, in general, *not* checked for unreasonable
</span><span class="comment">*</span><span class="comment">  values.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Arguments
</span><span class="comment">*</span><span class="comment">  =========
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  IJOB    (input) INTEGER
</span><span class="comment">*</span><span class="comment">          Specifies what is to be done:
</span><span class="comment">*</span><span class="comment">          = 1:  Compute NAB for the initial intervals.
</span><span class="comment">*</span><span class="comment">          = 2:  Perform bisection iteration to find eigenvalues of T.
</span><span class="comment">*</span><span class="comment">          = 3:  Perform bisection iteration to invert N(w), i.e.,
</span><span class="comment">*</span><span class="comment">                to find a point which has a specified number of
</span><span class="comment">*</span><span class="comment">                eigenvalues of T to its left.
</span><span class="comment">*</span><span class="comment">          Other values will cause <a name="SLAEBZ.70"></a><a href="slaebz.f.html#SLAEBZ.1">SLAEBZ</a> to return with INFO=-1.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  NITMAX  (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The maximum number of &quot;levels&quot; of bisection to be
</span><span class="comment">*</span><span class="comment">          performed, i.e., an interval of width W will not be made
</span><span class="comment">*</span><span class="comment">          smaller than 2^(-NITMAX) * W.  If not all intervals
</span><span class="comment">*</span><span class="comment">          have converged after NITMAX iterations, then INFO is set
</span><span class="comment">*</span><span class="comment">          to the number of non-converged intervals.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  N       (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The dimension n of the tridiagonal matrix T.  It must be at
</span><span class="comment">*</span><span class="comment">          least 1.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  MMAX    (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The maximum number of intervals.  If more than MMAX intervals
</span><span class="comment">*</span><span class="comment">          are generated, then <a name="SLAEBZ.85"></a><a href="slaebz.f.html#SLAEBZ.1">SLAEBZ</a> will quit with INFO=MMAX+1.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  MINP    (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The initial number of intervals.  It may not be greater than
</span><span class="comment">*</span><span class="comment">          MMAX.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  NBMIN   (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The smallest number of intervals that should be processed
</span><span class="comment">*</span><span class="comment">          using a vector loop.  If zero, then only the scalar loop
</span><span class="comment">*</span><span class="comment">          will be used.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  ABSTOL  (input) REAL
</span><span class="comment">*</span><span class="comment">          The minimum (absolute) width of an interval.  When an
</span><span class="comment">*</span><span class="comment">          interval is narrower than ABSTOL, or than RELTOL times the
</span><span class="comment">*</span><span class="comment">          larger (in magnitude) endpoint, then it is considered to be
</span><span class="comment">*</span><span class="comment">          sufficiently small, i.e., converged.  This must be at least
</span><span class="comment">*</span><span class="comment">          zero.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  RELTOL  (input) REAL
</span><span class="comment">*</span><span class="comment">          The minimum relative width of an interval.  When an interval
</span><span class="comment">*</span><span class="comment">          is narrower than ABSTOL, or than RELTOL times the larger (in
</span><span class="comment">*</span><span class="comment">          magnitude) endpoint, then it is considered to be
</span><span class="comment">*</span><span class="comment">          sufficiently small, i.e., converged.  Note: this should
</span><span class="comment">*</span><span class="comment">          always be at least radix*machine epsilon.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  PIVMIN  (input) REAL
</span><span class="comment">*</span><span class="comment">          The minimum absolute value of a &quot;pivot&quot; in the Sturm
</span><span class="comment">*</span><span class="comment">          sequence loop.  This *must* be at least  max |e(j)**2| *
</span><span class="comment">*</span><span class="comment">          safe_min  and at least safe_min, where safe_min is at least
</span><span class="comment">*</span><span class="comment">          the smallest number that can divide one without overflow.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  D       (input) REAL array, dimension (N)
</span><span class="comment">*</span><span class="comment">          The diagonal elements of the tridiagonal matrix T.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  E       (input) REAL array, dimension (N)
</span><span class="comment">*</span><span class="comment">          The offdiagonal elements of the tridiagonal matrix T in
</span><span class="comment">*</span><span class="comment">          positions 1 through N-1.  E(N) is arbitrary.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  E2      (input) REAL array, dimension (N)
</span><span class="comment">*</span><span class="comment">          The squares of the offdiagonal elements of the tridiagonal
</span><span class="comment">*</span><span class="comment">          matrix T.  E2(N) is ignored.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  NVAL    (input/output) INTEGER array, dimension (MINP)
</span><span class="comment">*</span><span class="comment">          If IJOB=1 or 2, not referenced.
</span><span class="comment">*</span><span class="comment">          If IJOB=3, the desired values of N(w).  The elements of NVAL
</span><span class="comment">*</span><span class="comment">          will be reordered to correspond with the intervals in AB.
</span><span class="comment">*</span><span class="comment">          Thus, NVAL(j) on output will not, in general be the same as
</span><span class="comment">*</span><span class="comment">          NVAL(j) on input, but it will correspond with the interval
</span><span class="comment">*</span><span class="comment">          (AB(j,1),AB(j,2)] on output.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  AB      (input/output) REAL array, dimension (MMAX,2)
</span><span class="comment">*</span><span class="comment">          The endpoints of the intervals.  AB(j,1) is  a(j), the left
</span><span class="comment">*</span><span class="comment">          endpoint of the j-th interval, and AB(j,2) is b(j), the
</span><span class="comment">*</span><span class="comment">          right endpoint of the j-th interval.  The input intervals
</span><span class="comment">*</span><span class="comment">          will, in general, be modified, split, and reordered by the
</span><span class="comment">*</span><span class="comment">          calculation.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  C       (input/output) REAL array, dimension (MMAX)
</span><span class="comment">*</span><span class="comment">          If IJOB=1, ignored.
</span><span class="comment">*</span><span class="comment">          If IJOB=2, workspace.
</span><span class="comment">*</span><span class="comment">          If IJOB=3, then on input C(j) should be initialized to the
</span><span class="comment">*</span><span class="comment">          first search point in the binary search.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  MOUT    (output) INTEGER
</span><span class="comment">*</span><span class="comment">          If IJOB=1, the number of eigenvalues in the intervals.
</span><span class="comment">*</span><span class="comment">          If IJOB=2 or 3, the number of intervals output.
</span><span class="comment">*</span><span class="comment">          If IJOB=3, MOUT will equal MINP.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  NAB     (input/output) INTEGER array, dimension (MMAX,2)
</span><span class="comment">*</span><span class="comment">          If IJOB=1, then on output NAB(i,j) will be set to N(AB(i,j)).
</span><span class="comment">*</span><span class="comment">          If IJOB=2, then on input, NAB(i,j) should be set.  It must
</span><span class="comment">*</span><span class="comment">             satisfy the condition:
</span><span class="comment">*</span><span class="comment">             N(AB(i,1)) &lt;= NAB(i,1) &lt;= NAB(i,2) &lt;= N(AB(i,2)),
</span><span class="comment">*</span><span class="comment">             which means that in interval i only eigenvalues
</span><span class="comment">*</span><span class="comment">             NAB(i,1)+1,...,NAB(i,2) will be considered.  Usually,
</span><span class="comment">*</span><span class="comment">             NAB(i,j)=N(AB(i,j)), from a previous call to <a name="SLAEBZ.160"></a><a href="slaebz.f.html#SLAEBZ.1">SLAEBZ</a> with
</span><span class="comment">*</span><span class="comment">             IJOB=1.
</span><span class="comment">*</span><span class="comment">             On output, NAB(i,j) will contain
</span><span class="comment">*</span><span class="comment">             max(na(k),min(nb(k),N(AB(i,j)))), where k is the index of
</span><span class="comment">*</span><span class="comment">             the input interval that the output interval
</span><span class="comment">*</span><span class="comment">             (AB(j,1),AB(j,2)] came from, and na(k) and nb(k) are the
</span><span class="comment">*</span><span class="comment">             the input values of NAB(k,1) and NAB(k,2).
</span><span class="comment">*</span><span class="comment">          If IJOB=3, then on output, NAB(i,j) contains N(AB(i,j)),
</span><span class="comment">*</span><span class="comment">             unless N(w) &gt; NVAL(i) for all search points  w , in which
</span><span class="comment">*</span><span class="comment">             case NAB(i,1) will not be modified, i.e., the output
</span><span class="comment">*</span><span class="comment">             value will be the same as the input value (modulo
</span><span class="comment">*</span><span class="comment">             reorderings -- see NVAL and AB), or unless N(w) &lt; NVAL(i)
</span><span class="comment">*</span><span class="comment">             for all search points  w , in which case NAB(i,2) will
</span><span class="comment">*</span><span class="comment">             not be modified.  Normally, NAB should be set to some
</span><span class="comment">*</span><span class="comment">             distinctive value(s) before <a name="SLAEBZ.174"></a><a href="slaebz.f.html#SLAEBZ.1">SLAEBZ</a> is called.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  WORK    (workspace) REAL array, dimension (MMAX)
</span><span class="comment">*</span><span class="comment">          Workspace.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  IWORK   (workspace) INTEGER array, dimension (MMAX)
</span><span class="comment">*</span><span class="comment">          Workspace.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  INFO    (output) INTEGER
</span><span class="comment">*</span><span class="comment">          = 0:       All intervals converged.
</span><span class="comment">*</span><span class="comment">          = 1--MMAX: The last INFO intervals did not converge.
</span><span class="comment">*</span><span class="comment">          = MMAX+1:  More than MMAX intervals were generated.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Further Details
</span><span class="comment">*</span><span class="comment">  ===============
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">      This routine is intended to be called only by other LAPACK
</span><span class="comment">*</span><span class="comment">  routines, thus the interface is less user-friendly.  It is intended
</span><span class="comment">*</span><span class="comment">  for two purposes:
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  (a) finding eigenvalues.  In this case, <a name="SLAEBZ.194"></a><a href="slaebz.f.html#SLAEBZ.1">SLAEBZ</a> should have one or
</span><span class="comment">*</span><span class="comment">      more initial intervals set up in AB, and <a name="SLAEBZ.195"></a><a href="slaebz.f.html#SLAEBZ.1">SLAEBZ</a> should be called
</span><span class="comment">*</span><span class="comment">      with IJOB=1.  This sets up NAB, and also counts the eigenvalues.
</span><span class="comment">*</span><span class="comment">      Intervals with no eigenvalues would usually be thrown out at
</span><span class="comment">*</span><span class="comment">      this point.  Also, if not all the eigenvalues in an interval i
</span><span class="comment">*</span><span class="comment">      are desired, NAB(i,1) can be increased or NAB(i,2) decreased.
</span><span class="comment">*</span><span class="comment">      For example, set NAB(i,1)=NAB(i,2)-1 to get the largest
</span><span class="comment">*</span><span class="comment">      eigenvalue.  <a name="SLAEBZ.201"></a><a href="slaebz.f.html#SLAEBZ.1">SLAEBZ</a> is then called with IJOB=2 and MMAX
</span><span class="comment">*</span><span class="comment">      no smaller than the value of MOUT returned by the call with
</span><span class="comment">*</span><span class="comment">      IJOB=1.  After this (IJOB=2) call, eigenvalues NAB(i,1)+1
</span><span class="comment">*</span><span class="comment">      through NAB(i,2) are approximately AB(i,1) (or AB(i,2)) to the
</span><span class="comment">*</span><span class="comment">      tolerance specified by ABSTOL and RELTOL.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  (b) finding an interval (a',b'] containing eigenvalues w(f),...,w(l).
</span><span class="comment">*</span><span class="comment">      In this case, start with a Gershgorin interval  (a,b).  Set up
</span><span class="comment">*</span><span class="comment">      AB to contain 2 search intervals, both initially (a,b).  One
</span><span class="comment">*</span><span class="comment">      NVAL element should contain  f-1  and the other should contain  l
</span><span class="comment">*</span><span class="comment">      , while C should contain a and b, resp.  NAB(i,1) should be -1
</span><span class="comment">*</span><span class="comment">      and NAB(i,2) should be N+1, to flag an error if the desired
</span><span class="comment">*</span><span class="comment">      interval does not lie in (a,b).  <a name="SLAEBZ.213"></a><a href="slaebz.f.html#SLAEBZ.1">SLAEBZ</a> is then called with
</span><span class="comment">*</span><span class="comment">      IJOB=3.  On exit, if w(f-1) &lt; w(f), then one of the intervals --
</span><span class="comment">*</span><span class="comment">      j -- will have AB(j,1)=AB(j,2) and NAB(j,1)=NAB(j,2)=f-1, while
</span><span class="comment">*</span><span class="comment">      if, to the specified tolerance, w(f-k)=...=w(f+r), k &gt; 0 and r
</span><span class="comment">*</span><span class="comment">      &gt;= 0, then the interval will have  N(AB(j,1))=NAB(j,1)=f-k and
</span><span class="comment">*</span><span class="comment">      N(AB(j,2))=NAB(j,2)=f+r.  The cases w(l) &lt; w(l+1) and
</span><span class="comment">*</span><span class="comment">      w(l-r)=...=w(l+k) are handled similarly.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  =====================================================================
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Parameters ..
</span>      REAL               ZERO, TWO, HALF
      PARAMETER          ( ZERO = 0.0E0, TWO = 2.0E0,
     $                   HALF = 1.0E0 / TWO )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      INTEGER            ITMP1, ITMP2, J, JI, JIT, JP, KF, KFNEW, KL,
     $                   KLNEW
      REAL               TMP1, TMP2
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          ABS, MAX, MIN
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Executable Statements ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Check for Errors
</span><span class="comment">*</span><span class="comment">
</span>      INFO = 0
      IF( IJOB.LT.1 .OR. IJOB.GT.3 ) THEN
         INFO = -1
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Initialize NAB
</span><span class="comment">*</span><span class="comment">
</span>      IF( IJOB.EQ.1 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Compute the number of eigenvalues in the initial intervals.
</span><span class="comment">*</span><span class="comment">
</span>         MOUT = 0
CDIR$ NOVECTOR
         DO 30 JI = 1, MINP
            DO 20 JP = 1, 2
               TMP1 = D( 1 ) - AB( JI, JP )
               IF( ABS( TMP1 ).LT.PIVMIN )
     $            TMP1 = -PIVMIN
               NAB( JI, JP ) = 0
               IF( TMP1.LE.ZERO )
     $            NAB( JI, JP ) = 1
<span class="comment">*</span><span class="comment">
</span>               DO 10 J = 2, N
                  TMP1 = D( J ) - E2( J-1 ) / TMP1 - AB( JI, JP )
                  IF( ABS( TMP1 ).LT.PIVMIN )
     $               TMP1 = -PIVMIN
                  IF( TMP1.LE.ZERO )
     $               NAB( JI, JP ) = NAB( JI, JP ) + 1
   10          CONTINUE
   20       CONTINUE
            MOUT = MOUT + NAB( JI, 2 ) - NAB( JI, 1 )
   30    CONTINUE
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Initialize for loop
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     KF and KL have the following meaning:
</span><span class="comment">*</span><span class="comment">        Intervals 1,...,KF-1 have converged.
</span><span class="comment">*</span><span class="comment">        Intervals KF,...,KL  still need to be refined.
</span><span class="comment">*</span><span class="comment">
</span>      KF = 1
      KL = MINP
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     If IJOB=2, initialize C.
</span><span class="comment">*</span><span class="comment">     If IJOB=3, use the user-supplied starting point.
</span><span class="comment">*</span><span class="comment">
</span>      IF( IJOB.EQ.2 ) THEN
         DO 40 JI = 1, MINP
            C( JI ) = HALF*( AB( JI, 1 )+AB( JI, 2 ) )
   40    CONTINUE
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Iteration loop
</span><span class="comment">*</span><span class="comment">
</span>      DO 130 JIT = 1, NITMAX
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Loop over intervals
</span><span class="comment">*</span><span class="comment">
</span>         IF( KL-KF+1.GE.NBMIN .AND. NBMIN.GT.0 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Begin of Parallel Version of the loop
</span><span class="comment">*</span><span class="comment">
</span>            DO 60 JI = KF, KL
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Compute N(c), the number of eigenvalues less than c
</span><span class="comment">*</span><span class="comment">
</span>               WORK( JI ) = D( 1 ) - C( JI )
               IWORK( JI ) = 0
               IF( WORK( JI ).LE.PIVMIN ) THEN
                  IWORK( JI ) = 1
                  WORK( JI ) = MIN( WORK( JI ), -PIVMIN )
               END IF
<span class="comment">*</span><span class="comment">
</span>               DO 50 J = 2, N
                  WORK( JI ) = D( J ) - E2( J-1 ) / WORK( JI ) - C( JI )
                  IF( WORK( JI ).LE.PIVMIN ) THEN
                     IWORK( JI ) = IWORK( JI ) + 1
                     WORK( JI ) = MIN( WORK( JI ), -PIVMIN )
                  END IF
   50          CONTINUE
   60       CONTINUE
<span class="comment">*</span><span class="comment">
</span>            IF( IJOB.LE.2 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              IJOB=2: Choose all intervals containing eigenvalues.
</span><span class="comment">*</span><span class="comment">
</span>               KLNEW = KL
               DO 70 JI = KF, KL
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Insure that N(w) is monotone
</span><span class="comment">*</span><span class="comment">
</span>                  IWORK( JI ) = MIN( NAB( JI, 2 ),
     $                          MAX( NAB( JI, 1 ), IWORK( JI ) ) )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Update the Queue -- add intervals if both halves
</span><span class="comment">*</span><span class="comment">                 contain eigenvalues.
</span><span class="comment">*</span><span class="comment">
</span>                  IF( IWORK( JI ).EQ.NAB( JI, 2 ) ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    No eigenvalue in the upper interval:
</span><span class="comment">*</span><span class="comment">                    just use the lower interval.
</span><span class="comment">*</span><span class="comment">
</span>                     AB( JI, 2 ) = C( JI )
<span class="comment">*</span><span class="comment">
</span>                  ELSE IF( IWORK( JI ).EQ.NAB( JI, 1 ) ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    No eigenvalue in the lower interval:
</span><span class="comment">*</span><span class="comment">                    just use the upper interval.
</span><span class="comment">*</span><span class="comment">
</span>                     AB( JI, 1 ) = C( JI )
                  ELSE
                     KLNEW = KLNEW + 1
                     IF( KLNEW.LE.MMAX ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                       Eigenvalue in both intervals -- add upper to
</span><span class="comment">*</span><span class="comment">                       queue.
</span><span class="comment">*</span><span class="comment">
</span>                        AB( KLNEW, 2 ) = AB( JI, 2 )
                        NAB( KLNEW, 2 ) = NAB( JI, 2 )
                        AB( KLNEW, 1 ) = C( JI )
                        NAB( KLNEW, 1 ) = IWORK( JI )
                        AB( JI, 2 ) = C( JI )
                        NAB( JI, 2 ) = IWORK( JI )
                     ELSE
                        INFO = MMAX + 1
                     END IF
                  END IF
   70          CONTINUE
               IF( INFO.NE.0 )
     $            RETURN
               KL = KLNEW
            ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              IJOB=3: Binary search.  Keep only the interval containing
</span><span class="comment">*</span><span class="comment">                      w   s.t. N(w) = NVAL
</span><span class="comment">*</span><span class="comment">
</span>               DO 80 JI = KF, KL
                  IF( IWORK( JI ).LE.NVAL( JI ) ) THEN
                     AB( JI, 1 ) = C( JI )
                     NAB( JI, 1 ) = IWORK( JI )
                  END IF
                  IF( IWORK( JI ).GE.NVAL( JI ) ) THEN
                     AB( JI, 2 ) = C( JI )
                     NAB( JI, 2 ) = IWORK( JI )
                  END IF
   80          CONTINUE
            END IF
<span class="comment">*</span><span class="comment">
</span>         ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           End of Parallel Version of the loop
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Begin of Serial Version of the loop
</span><span class="comment">*</span><span class="comment">
</span>            KLNEW = KL
            DO 100 JI = KF, KL
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Compute N(w), the number of eigenvalues less than w
</span><span class="comment">*</span><span class="comment">
</span>               TMP1 = C( JI )
               TMP2 = D( 1 ) - TMP1
               ITMP1 = 0
               IF( TMP2.LE.PIVMIN ) THEN
                  ITMP1 = 1
                  TMP2 = MIN( TMP2, -PIVMIN )
               END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              A series of compiler directives to defeat vectorization
</span><span class="comment">*</span><span class="comment">              for the next loop
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">$PL$ CMCHAR=' '
</span>CDIR$          NEXTSCALAR
C$DIR          SCALAR
CDIR$          NEXT SCALAR
CVD$L          NOVECTOR
CDEC$          NOVECTOR
CVD$           NOVECTOR
<span class="comment">*</span><span class="comment">VDIR          NOVECTOR
</span><span class="comment">*</span><span class="comment">VOCL          LOOP,SCALAR
</span>CIBM           PREFER SCALAR
<span class="comment">*</span><span class="comment">$PL$ CMCHAR='*'
</span><span class="comment">*</span><span class="comment">
</span>               DO 90 J = 2, N
                  TMP2 = D( J ) - E2( J-1 ) / TMP2 - TMP1
                  IF( TMP2.LE.PIVMIN ) THEN
                     ITMP1 = ITMP1 + 1
                     TMP2 = MIN( TMP2, -PIVMIN )
                  END IF
   90          CONTINUE
<span class="comment">*</span><span class="comment">
</span>               IF( IJOB.LE.2 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 IJOB=2: Choose all intervals containing eigenvalues.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Insure that N(w) is monotone
</span><span class="comment">*</span><span class="comment">
</span>                  ITMP1 = MIN( NAB( JI, 2 ),
     $                    MAX( NAB( JI, 1 ), ITMP1 ) )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 Update the Queue -- add intervals if both halves
</span><span class="comment">*</span><span class="comment">                 contain eigenvalues.
</span><span class="comment">*</span><span class="comment">
</span>                  IF( ITMP1.EQ.NAB( JI, 2 ) ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    No eigenvalue in the upper interval:
</span><span class="comment">*</span><span class="comment">                    just use the lower interval.
</span><span class="comment">*</span><span class="comment">
</span>                     AB( JI, 2 ) = TMP1
<span class="comment">*</span><span class="comment">
</span>                  ELSE IF( ITMP1.EQ.NAB( JI, 1 ) ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    No eigenvalue in the lower interval:
</span><span class="comment">*</span><span class="comment">                    just use the upper interval.
</span><span class="comment">*</span><span class="comment">
</span>                     AB( JI, 1 ) = TMP1
                  ELSE IF( KLNEW.LT.MMAX ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                    Eigenvalue in both intervals -- add upper to queue.
</span><span class="comment">*</span><span class="comment">
</span>                     KLNEW = KLNEW + 1
                     AB( KLNEW, 2 ) = AB( JI, 2 )
                     NAB( KLNEW, 2 ) = NAB( JI, 2 )
                     AB( KLNEW, 1 ) = TMP1
                     NAB( KLNEW, 1 ) = ITMP1
                     AB( JI, 2 ) = TMP1
                     NAB( JI, 2 ) = ITMP1
                  ELSE
                     INFO = MMAX + 1
                     RETURN
                  END IF
               ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">                 IJOB=3: Binary search.  Keep only the interval
</span><span class="comment">*</span><span class="comment">                         containing  w  s.t. N(w) = NVAL
</span><span class="comment">*</span><span class="comment">
</span>                  IF( ITMP1.LE.NVAL( JI ) ) THEN
                     AB( JI, 1 ) = TMP1
                     NAB( JI, 1 ) = ITMP1
                  END IF
                  IF( ITMP1.GE.NVAL( JI ) ) THEN
                     AB( JI, 2 ) = TMP1
                     NAB( JI, 2 ) = ITMP1
                  END IF
               END IF
  100       CONTINUE
            KL = KLNEW
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           End of Serial Version of the loop
</span><span class="comment">*</span><span class="comment">
</span>         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Check for convergence
</span><span class="comment">*</span><span class="comment">
</span>         KFNEW = KF
         DO 110 JI = KF, KL
            TMP1 = ABS( AB( JI, 2 )-AB( JI, 1 ) )
            TMP2 = MAX( ABS( AB( JI, 2 ) ), ABS( AB( JI, 1 ) ) )
            IF( TMP1.LT.MAX( ABSTOL, PIVMIN, RELTOL*TMP2 ) .OR.
     $          NAB( JI, 1 ).GE.NAB( JI, 2 ) ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Converged -- Swap with position KFNEW,
</span><span class="comment">*</span><span class="comment">                           then increment KFNEW
</span><span class="comment">*</span><span class="comment">
</span>               IF( JI.GT.KFNEW ) THEN
                  TMP1 = AB( JI, 1 )
                  TMP2 = AB( JI, 2 )
                  ITMP1 = NAB( JI, 1 )
                  ITMP2 = NAB( JI, 2 )
                  AB( JI, 1 ) = AB( KFNEW, 1 )
                  AB( JI, 2 ) = AB( KFNEW, 2 )
                  NAB( JI, 1 ) = NAB( KFNEW, 1 )
                  NAB( JI, 2 ) = NAB( KFNEW, 2 )
                  AB( KFNEW, 1 ) = TMP1
                  AB( KFNEW, 2 ) = TMP2
                  NAB( KFNEW, 1 ) = ITMP1
                  NAB( KFNEW, 2 ) = ITMP2
                  IF( IJOB.EQ.3 ) THEN
                     ITMP1 = NVAL( JI )
                     NVAL( JI ) = NVAL( KFNEW )
                     NVAL( KFNEW ) = ITMP1
                  END IF
               END IF
               KFNEW = KFNEW + 1
            END IF
  110    CONTINUE
         KF = KFNEW
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Choose Midpoints
</span><span class="comment">*</span><span class="comment">
</span>         DO 120 JI = KF, KL
            C( JI ) = HALF*( AB( JI, 1 )+AB( JI, 2 ) )
  120    CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        If no more intervals to refine, quit.
</span><span class="comment">*</span><span class="comment">
</span>         IF( KF.GT.KL )
     $      GO TO 140
  130 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Converged
</span><span class="comment">*</span><span class="comment">
</span>  140 CONTINUE
      INFO = MAX( KL+1-KF, 0 )
      MOUT = KL
<span class="comment">*</span><span class="comment">
</span>      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="SLAEBZ.549"></a><a href="slaebz.f.html#SLAEBZ.1">SLAEBZ</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>
